#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "MSTK.h"


#ifdef __cplusplus
extern "C" {
#endif

  /* 
     This function adds ghost entities and label PUTILIZED entities

     for face(surface mesh) or region(volume mesh), add ghost and 
     label utilized, then add and label corresponding vertices.

     for boundary vertices, label utilized and ghost.

     ring: ghost cell size, 
     0: only ghost vertex on processor boundary
     1: 1-ring of ghost face(surface mesh) or region(volume mesh)

     must call PMESH_BuildPBoundary() first
          
  */
int MESH_Vol_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int part_no, int ring);
int MESH_Surf_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int part_no, int ring);


int MESH_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int part_no, int ring) {
  int nf, nr, idx;
  MVertex_ptr local_mv;

  /* if 0 ring, only mark processor boundary vertex */
  if (ring == 0) {
    idx = 0;
    while((local_mv = MESH_Next_Vertex(submesh,&idx))) {
      if(MV_PType(local_mv) == PBOUNDARY) {
	if(MV_MasterParID(local_mv) == part_no)
	  MV_Set_PType(local_mv,PUTILIZED);
	else
	  MV_Set_PType(local_mv,PGHOST);
      }
    }
    return 1;
  }

  /* basic mesh information */
  nf = MESH_Num_Faces(submesh);
  nr = MESH_Num_Regions(submesh);

  if (nr)
    MESH_Vol_AddGhost(mesh,submesh,part_no,ring);
  else if(nf) 
    MESH_Surf_AddGhost(mesh,submesh,part_no,ring);
  else {
    MSTK_Report("MESH_AddGhost()","This is not a valid mstk file",ERROR);
    exit(-1);
  }
  return 1;
}

int MESH_Surf_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int part_no, int ring) {
  int iv, i, j, idx, vertex_id, face_id, mkvid, mkfid;
  int nv, nf, nvf, nfv;
  MVertex_ptr mv, mv_new, *fverts, global_mv, local_mv;
  MFace_ptr local_mf, global_mf, mf_new;
  List_ptr global_mvfaces, global_mfvertices,local_mfvertices;
  List_ptr par_vertices, par_faces;
  int *v_g2l, *f_g2l;
  /* basic mesh information */
  nv = MESH_Num_Vertices(submesh);
  nf = MESH_Num_Faces(submesh);

  /* store the map between the vertex global ID to ID in submesh */
  v_g2l = (int *)MSTK_malloc((MESH_Num_Vertices(mesh)+1)*sizeof(int));
  /* store the map between the face global ID to ID in submesh */
  f_g2l = (int *)MSTK_malloc((MESH_Num_Faces(mesh)+1)*sizeof(int));

  /* Mark the list of global entities in this submesh */
  par_vertices = List_New(nv);
  mkvid = MSTK_GetMarker();
  idx = 0;
  while((local_mv = MESH_Next_Vertex(submesh,&idx))) {
    vertex_id = MV_GlobalID(local_mv);
    /* fill the map */
    v_g2l[vertex_id] = MV_ID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    List_Add(par_vertices,global_mv);
    MEnt_Mark(global_mv,mkvid);
  }


  par_faces = List_New(nf);
  mkfid = MSTK_GetMarker();  
  idx = 0;
  while((local_mf = MESH_Next_Face(submesh,&idx))) {
    face_id = MF_GlobalID(local_mf);
    /* fill the map */
    f_g2l[face_id] = MF_ID(local_mf);
    global_mf = MESH_FaceFromID(mesh,face_id);
    List_Add(par_faces,global_mf);
    MEnt_Mark(global_mf,mkfid);
  }

  /*Loop through sub vertices, and see its neighboring faces */
  for(iv = 0; iv < nv; iv++) {
    /* get global face id and loop through global mesh
     here, can not use List par_vertices, it will be dynamically added */
    local_mv = MESH_Vertex(submesh,iv);
    /* only check vertices on the processor boundary */
    if(MV_PType(local_mv) != PBOUNDARY)
      continue;
    vertex_id = MV_GlobalID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    /* get global neighbor faces of this vertex */
    global_mvfaces = MV_Faces(global_mv);
    nvf = List_Num_Entries(global_mvfaces);
    for(i = 0; i < nvf; i++) {
      global_mf = List_Entry(global_mvfaces,i);
      /* if the face is already marked */
      if (MEnt_IsMarked(global_mf,mkfid)) {
	local_mf = MESH_FaceFromID(submesh,f_g2l[MF_GlobalID(global_mf)]);
	/* mark PUTILIZED face */
	if(MF_PType(local_mf) == PINTERIOR) {
	  MF_Set_PType(local_mf,PUTILIZED);
	  local_mfvertices = MF_Vertices(local_mf,1,0);
	  for(j = 0; j < List_Num_Entries(local_mfvertices); j++) {
	    mv = List_Entry(local_mfvertices,j);
	    /* mark PUTILIZED vertex */
	    if( (MV_PType(mv) != PBOUNDARY) && (MV_MasterParID(mv) == part_no) )
	      MV_Set_PType(mv,PUTILIZED);
	  }
	  List_Delete(local_mfvertices);
	}
	continue;
      }
      /* if it is not marked, add this face into submesh */
      global_mfvertices = MF_Vertices(global_mf,1,0);
      nfv = List_Num_Entries(global_mfvertices);
      fverts = (MVertex_ptr *)MSTK_malloc((nfv)*sizeof(MVertex_ptr));
      /* loop through the global vertices of this face */
      for(j = 0; j < nfv; j++) {
	global_mv = List_Entry(global_mfvertices,j);
	/* if this vertex is already marked, use its local counterpart */
	if(MEnt_IsMarked(global_mv,mkvid)) {
	  mv = MESH_VertexFromID(submesh,v_g2l[MV_GlobalID(global_mv)]);
	  fverts[j] = mv;
	}
	/* otherwise, add this vertex into submesh */
	if(!MEnt_IsMarked(global_mv,mkvid)) {
	  List_Add(par_vertices,global_mv);
	  MEnt_Mark(global_mv,mkvid);
	  mv_new = MV_New(submesh);
	  PMV_Duplicate(global_mv,mv_new);
	  /* ptype is not copied in PMV_Copy */
	  MV_Set_PType(mv_new,PGHOST);
	  fverts[j] = mv_new;
	  /* fill the map */
	  v_g2l[MV_GlobalID(global_mv)] = MV_ID(mv_new);
	}
      }
      /* add face */
      List_Add(par_faces,global_mf);
      MEnt_Mark(global_mf,mkfid);
      mf_new = MF_New(submesh);
      MF_Set_Vertices(mf_new,nfv,fverts);
      MF_Set_PType(mf_new,PGHOST);
      PMF_Duplicate(global_mf,mf_new);
      /* fill the map */
      f_g2l[MF_GlobalID(global_mf)] = MF_ID(mf_new);

      if (fverts) MSTK_free(fverts);
      List_Delete(global_mfvertices);
    }
    List_Delete(global_mvfaces);
    /* and mark the PBOUNDARY as PUTILIZED or PGHOST */
    if(MV_MasterParID(local_mv) == part_no)
      MV_Set_PType(local_mv,PUTILIZED);
    else
      MV_Set_PType(local_mv,PGHOST);
  }

  List_Unmark(par_vertices,mkvid);
  List_Delete(par_vertices);

  List_Unmark(par_faces,mkfid);
  List_Delete(par_faces);

  MSTK_FreeMarker(mkvid);
  MSTK_FreeMarker(mkfid);

  MSTK_free(v_g2l);  
  MSTK_free(f_g2l);  


  return 1;
}

int MESH_Vol_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int part_no, int ring) {
  int iv, i, j, idx, vertex_id, region_id, mkvid, mkrid;
  int nv, nr, nvr, nrv;
  MVertex_ptr mv, mv_new, *rverts, global_mv, local_mv;
  MRegion_ptr local_mr, global_mr, mr_new;
  List_ptr global_mvregions, global_mrvertices,local_mrvertices;
  List_ptr par_vertices, par_regions;
  int *v_g2l, *r_g2l;
  /* basic mesh information */
  nv = MESH_Num_Vertices(submesh);
  nr = MESH_Num_Regions(submesh);

  /* store the map between the vertex global ID to ID in submesh */
  v_g2l = (int *)MSTK_malloc((MESH_Num_Vertices(mesh)+1)*sizeof(int));
  /* store the map between the region global ID to ID in submesh */
  r_g2l = (int *)MSTK_malloc((MESH_Num_Regions(mesh)+1)*sizeof(int));

  /* Mark the list of global entities in this submesh */
  par_vertices = List_New(nv);
  mkvid = MSTK_GetMarker();
  idx = 0;
  while((local_mv = MESH_Next_Vertex(submesh,&idx))) {
    vertex_id = MV_GlobalID(local_mv);
    /* fill the map */
    v_g2l[vertex_id] = MV_ID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    List_Add(par_vertices,global_mv);
    MEnt_Mark(global_mv,mkvid);
  }


  par_regions = List_New(nr);
  mkrid = MSTK_GetMarker();  
  idx = 0;
  while((local_mr = MESH_Next_Region(submesh,&idx))) {
    region_id = MR_GlobalID(local_mr);
    /* fill the map */
    r_g2l[region_id] = MR_ID(local_mr);
    global_mr = MESH_RegionFromID(mesh,region_id);
    List_Add(par_regions,global_mr);
    MEnt_Mark(global_mr,mkrid);
  }

  /*Loop through sub vertices, and see its neighboring regions */
  for(iv = 0; iv < nv; iv++) {
    /* get global region id and loop through global mesh
     here, can not use List par_vertices, it will be dynamically added */
    local_mv = MESH_Vertex(submesh,iv);
    /* only check vertices on the processor boundary */
    if(MV_PType(local_mv) != PBOUNDARY)
      continue;
    vertex_id = MV_GlobalID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    /* get global neighbor regions of this vertex */
    global_mvregions = MV_Regions(global_mv);
    nvr = List_Num_Entries(global_mvregions);
    for(i = 0; i < nvr; i++) {
      global_mr = List_Entry(global_mvregions,i);
      /* if the region is already marked */
      if (MEnt_IsMarked(global_mr,mkrid)) {
	local_mr = MESH_RegionFromID(submesh,r_g2l[MR_GlobalID(global_mr)]);
	/* mark PUTILIZED region */
	if(MR_PType(local_mr) == PINTERIOR) {
	  MR_Set_PType(local_mr,PUTILIZED);
	  local_mrvertices = MR_Vertices(local_mr);
	  for(j = 0; j < List_Num_Entries(local_mrvertices); j++) {
	    mv = List_Entry(local_mrvertices,j);
	    /* mark PUTILIZED vertex */
	    if( (MV_PType(mv) != PBOUNDARY) && (MV_MasterParID(mv) == part_no) )
	      MV_Set_PType(mv,PUTILIZED);
	  }
	  List_Delete(local_mrvertices);
	}
	continue;
      }
      /* if it is not marked, add this region into submesh */
      global_mrvertices = MR_Vertices(global_mr);
      nrv = List_Num_Entries(global_mrvertices);
      rverts = (MVertex_ptr *)MSTK_malloc((nrv)*sizeof(MVertex_ptr));
      /* loop through the global vertices of this region */
      for(j = 0; j < nrv; j++) {
	global_mv = List_Entry(global_mrvertices,j);
	/* if this vertex is already marked, use its local counterpart */
	if(MEnt_IsMarked(global_mv,mkvid)) {
	  mv = MESH_VertexFromID(submesh,v_g2l[MV_GlobalID(global_mv)]);
	  rverts[j] = mv;
	}
	/* otherwise, add this vertex into submesh */
	if(!MEnt_IsMarked(global_mv,mkvid)) {
	  List_Add(par_vertices,global_mv);
	  MEnt_Mark(global_mv,mkvid);
	  mv_new = MV_New(submesh);
	  PMV_Duplicate(global_mv,mv_new);
	  /* ptype is not copied in PMV_Copy */
	  MV_Set_PType(mv_new,PGHOST);
	  rverts[j] = mv_new;
	  /* fill the map */
	  v_g2l[MV_GlobalID(global_mv)] = MV_ID(mv_new);
	}
      }
      /* add region */
      List_Add(par_regions,global_mr);
      MEnt_Mark(global_mr,mkrid);
      mr_new = MR_New(submesh);
      MR_Set_Vertices(mr_new,nrv,rverts,0,NULL);
      MR_Set_PType(mr_new,PGHOST);
      PMR_Duplicate(global_mr,mr_new);
      /* fill the map */
      r_g2l[MR_GlobalID(global_mr)] = MR_ID(mr_new);

      if (rverts) MSTK_free(rverts);
      List_Delete(global_mrvertices);
    }
    List_Delete(global_mvregions);
    /* and mark the PBOUNDARY as PUTILIZED or PGHOST */
    if(MV_MasterParID(local_mv) == part_no)
      MV_Set_PType(local_mv,PUTILIZED);
    else
      MV_Set_PType(local_mv,PGHOST);
  }

  List_Unmark(par_vertices,mkvid);
  List_Delete(par_vertices);

  List_Unmark(par_regions,mkrid);
  List_Delete(par_regions);

  MSTK_FreeMarker(mkvid);
  MSTK_FreeMarker(mkrid);

  MSTK_free(v_g2l);  
  MSTK_free(r_g2l);  


  return 1;
}
  
#ifdef __cplusplus
}
#endif

